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SUMMARY 

at the University of Texas at El Paso’s (UTEP’s) research solar pond. 


INTRODUCTION 

to , non-convecting salinity gradient solar pond, natural co nvection it 

: sr; 

SilJSSSStoZi to ensure tat <h« density 

of the fluid is continuously increasing. 

H* stability of the solar pond is a function of the thennal and solutal differ rat« 
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temperature and salinity profiles which can be expected over an operating period when estima g 

daily ambient temperatures and insolation values are provided. 
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Background 


Solar energy incident upon the surface of a shallow pond will be partially transmitted to the bottom of 
the pond where much of it is absorbed. This normally leads to natural convection as absorbed energy 
in die lower regions of the pond causes the fluid to thermally expand. Consequently, the thermal energy 
will be transported to the surface and released. In a non-convecting salinity gradient solar pond this 
convective process is suppressed by establishing an internal region with a salinity gradient of sufficient 
strength to ensure that the localized decreases in density due to thermal expansion are counter-balanced 
by die localized increases in density due to the higher salt concentrations. Thermal conduction will then 
become the dominant heat transfer release mechanism through die region of the pond where convection 
is suppressed. Due to the relatively low value of thermal conductivity for water (approximately 0.65 
W/m-°C) an "insulating'' layer develops which allows temperatures approaching 100 °C in the pond’s 
thermal storage zone. 

The typical solar pond is characterized by three regions (Figure 1); a convective surface region referred 
to as the Upper Convecting Zone (UCZ), a non-convective primary gradient region referred to as the 
Non-Convecting Zone (NCZ), and a convective thermal storage region referred to as the Lower 
Convecting Zone (LCZ). The thickness of each of these zones depends upon die season of the year and 
the manner in which the pond is being operated. For purposes of this analysis average values were 
derived from the research pond currently in operation at the University of Texas at El Paso (UTEP). The 
dimensions of this pond are 52.4 m across by 64.0 m long comprising a total surface area of 3,355 m : 
(about 5/6 of an acre). The walls of the solar pond slope inward at an approximate angle of 30° from 
horizontal. The depth of the pond is about 3.5 meters resulting in an approximate floor area of 2,500 
m 2 . 
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Figure 1. Schematic of a solar pond. 
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MODELING APPROACH 


period of time. 


nerivation of g puraW Conditions 

Emiier work wm performed by Hull showing fo« modeling *« hourly ' S 

STL.lts obuJdfor foe LCZ remperafore over an «teoded penod o f afody [1]. Sm* dmrmd 

dady t3E -? mrtjten? temperature were used for foe 

transient boundary conditions. 

A. interpolating function was derived for insolation on a horironul surf to f rom monthly avenges 
obtained over a 21 year period for El Paso, Texas 121. This is given as 


-6,3W*2,3W-sinplg^j 


(W/m 1 ) 


Eq. (1) 


where: 


n = day of the year. 


The portion of solar radiation incident on a horizontal surface that is reflected at Ae sur^ of Ae air- 
wJer Sace of *e pond is modeled by assuming that the incident radiation is d^ ^ r^on 
:^=rthe surface at a fixed angle of incidence. These assumptions were previously reviewed by 

fto\\ in a comparison of results from a detailed computer model with an J^^^CzS^^re 
assumption of direct beam radiation. The results showed that accurate estimates of the LCZ temp 
^T obteined when the fixed angle of incidence (i.e., A,) is calculated at solar noon 17 days before the 

autumnal equinox [11. 


The following equation is used to calculate the fixed angle of incidence used in this analysis [2]: 


a. - cos* 1 (sin/S *siny + cos/5 *0057 *cosw) 


Eq- (2) 


where: 


y 

fi 


u — 


latitude of El Paso, Texas (i.e., 31.5°). 

declination of sun on the 17th day preceding the autumnal equinox, 
hour angle at solar noon (i.e, t 0°). 


The angle of declination (/5) is the angle of the earth’s axis relative to the sun-earth line and is calculated 
for die 17th day preceding the autumnal equinox as follows. 


0 = 23.44° *sin 


(360 ‘(5-81)) 
365.25 


Eq. (3) 


where: 6 = the day of the year for September 6th (i.e., 249) 
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This gives a calculated value for a, of 25.66° which is the fixed angle of incidence used in this analysis. 
The associated fixed angle of refraction is calculated using Snell’s Law by the following relationship: 


p F * sin 


(sin(cr,)) 

n 


Eq. (4) 


Using 1.333 for the index of refraction, (n), for the pond’s salt water solution [3] Snell’s Law yields a 
value of 18.99° for foe fixed angle of refraction, p F . 

The surface reflectivity is then calculated using foe following relationship derived from Fresnel’s 
equations: 

. . .5 . fgfevgd - « ■ tm ’< <, r<y> £,.(5) 

* sin^a,.?,) tan’(a,+p f ) 


Applying foe calculated values for the fixed angle of incidence and refraction to this equation yields an 
effective surface reflectivity ( p *) of .021. 

The ambient temperature is calculated in a manner similar to that used in developing foe average daily 
insolation. Average daily ambient temperatures are obtained from a periodic interpolating function fitted 
to actual historical temperatures recorded over a thirty year period by foe National Oceanic and 
Atmospheric Administration. 


TJLn) 


63.55 + 19. 15 ‘Sin 


(2x(n - 105)) 
365.25 


(°C) 


Eq. (6) 


where: n = foe day of foe year. 

Pond Attenuation of Solar Radiation 

As foe radiation is transmitted through foe pond it is attenuated along it’s path by scattering and 
absorption. The shorter wavelengths of the solar radiation will be transmitted through foe pond with very 
little absorption, whereas the longer wavelengths (i.e., foe infrared portion of foe spectrum) will be 
absorbed within foe first few centimeters. Wavelengths longer than foe infrared wavelength (i.e., greater 
than Kftn) are reflected at foe pond’s surface. 

The absorption and scattering of foe solar radiation is generally represented by a transmission function 
from which foe radiation intensity at a given depth can be determined. One such function, developed by 
Hull, was derived from an extensive set of experimental absorption data measured for pure water. This 
data was reduced into a four part summation of exponential terms. 

t * E<*p ^ w 

where: r = foe fraction of solar radiation reaching a depth (x). 
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This reduced function correlates well with a more extensive 
data, never deviating by more than .01 in transmittance over 
used' for this transmission function are provided in Table I. 


forty part function, derived from the same 
a 2 meter interval. The parametric values 



Table 


Exponential Extinction Coefficients for Water 


0) 


Bj 


Hj (m’ 1 ) 


0.237 

0.193 

0.167 

0.179 


0.032 

0.450 

3.000 

35.000 


When compared with nKasurements taken at the UTEP research solas pondon April 20 k»W ^ HnU 
Four-Part Junction resulted in substantially higher transmittance vAies, which is to 

transmission Junction developed by Rabl-Nielson [4) have been presented in Figure 2 for comparison. 

Eq. (8) 




= - 066i 3 +.385x 2 -.768jc+.745 


^tahiiirv Criteria Conside raliQBS 

He NCZ has been subdivided into two sub-regions to recognize the t no« ‘ artagems. ^ 
is suppressed: 

( Pr*r ) .P.n B* W 

r (Pr+1) a 


where: 


G t = 
G s = 

& = 
a = 
Pr = 

r = 


the temperature gradient at various locations in the NCZ (°C/m), 

the salinity gradient at corresponding locations (%/m), 

the saline expansion coefficient for NaCl (m J /Kg), 

the thermal expansion coefficient (1/°C), . , 

Prandtl number (i.e., ratio of the kinematic viscosity to the thermal 

diffusivity), . „„ . . v rv n 

ratio of the saline diffusivity to the thermal diffusivity (».e., k s /k t ). 
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Figure 2. Comparative plot of transmission (unctions. 


Since the values for the parameters in the above equation are dependent on both temperature and salinity 
it is necessary to develop a functional relationship between the salinity gradient factor (i.e., SGF, and 
defined as SGF = [(Pr + r)/(Pr - l))*(/3/a) ) and temperature [6]. Using information published by the 
Department of Interior’s Office of Saline Water the following relationship was derived [7]: 

SGF(T ) = *~Ill 19 ?-?5 2 ) - (1. 702)7’ 48. 876 Eq. (10) 

7,812.5 


A plot of Equation 10 with the actual data for the SGF is presented in Figure 3. This relationship was 
then inserted into Equation 9 to obtain the necessary temperature and salinity gradient relationship that 
was used in establishing the internal threshold salinity gradient. 

The dynamic stability criterion that is employed at both the upper and lower NCZ interfaces is based on 
an empirical relationship derived from Neilsen and is referred to as the Neilsen boundary condition [8]. 
This empirical relationship is given as follows: 

G s -A-G? Eq.(ll) 

where: A is set at 28.0 (Kg/m 4 ) • (m/K) ® 

Development of Temperature Profile Equations 

The nodal equations assume the thermal exchanges that take place between the various zones within the 
solar pond can be adequately represented by a transient one-dimensional model. Thermal properties of 
the pond are assumed to be uniform within each of the pond’s zones. Using these assumptions, 
temperature changes, resulting from thermal exchanges within the pond, can then be represented by a 
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T#n*r«Tir* c DtorHS C 3 
■ Actual Data 


Figure 3. Comparative plot of the SGF interpolating function with actual data 


I^o^ionbi a system of one^dimensional foretarWifference ***** <«■“« * is postfve 
in the downward direction from the surface at x — 0.0 ). 

Since the UCZ is essentially uniform in temperature it is represented by a single node. It's 

tsi^^ST^mfog n^^S^^dien. L^t foe »«. 

JKSiSr^Sif fo d”fo of foe NCZ-LCZ interface, and beyond. is assumed fo be»mpl«ely 
ambient temperature exists. 


UCZ-Nodal Equation: 


Tyez - T,(f) 


Eq. (12) 


UCZ-NCZ Interface Nodal Equation: w-Fo 

Toj i * T 0 ; (1 -IFo-lFo • BO *2Fo • (T t J *Bi • T^) + — jj — EQ- (13) 
NCZ Interior Nodal Equations: fix *Fo 

T = T.j • (1 -2 Fo) * Fo • (T mj * T mj ) * — ^ — * Qi-xnj-h.mj’ Eq. (14) 
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NCZ-LCZ Interface Nodal Equation: 


2 Bi’Fo 


- T m y(l-2Fo- 2 Fo • Bi) ♦ 2 Fo • (T.. w *« • 7^, ) ♦ ^1 


(/. 




Eq. (15) 


LCZ Interior Nodal Equation: 

W, - r w .(i-2»).ii-(vv*2-W«) 


where: 


H 


dt'h c 

L ta m *' c e 


Eq. (16) 


LCZ-GRZ Interfacial Nodal Equation: 

V, = 7;, • (1 -2Fo-2Fo • B/) * 2Fo • (T^ * Bi • T^) Eq. (17) 


GRZ Interior Nodal Equation: 

T.J., ' T,y(\-2Fo) . Fo-U^.T^ Eq. (18) 


In the above equations and T |f7j are the bulk temperatures of the UCZ and LCZ at time (j), 
respectively, Fo is the Fourier Modulus, and Bi is the Biot modulus. The heat transfer coefficient is set 
at 263 W/m 2 -°C and is assumed to be the same at both interfaces. This value is derived from 
consideration of the empirical relationship developed for two horizontal plates. 


Development of Salinity Profile Equations 


The salinity profile is derived from the solution of a transient one dimensional partial differential equation 
wherein the rate of salinity change is presented as a function of the temperature dependent solutal 
diffusivity (K s ) and the rate of change of salinity with respect to depth in the gradient zone. After the 
salinity profile is calculated for each progression in temperature profile it is checked to ensure that die 
stability criteria at both the upper and lower boundaries as well as within the NCZ are satisfied. If the 
boundary stability criterion at either die upper or lower interfaces is not met the boundary is repositioned 
until the stability criterion is satisfied. If die internal stability criterion is not met the run is simply 
terminated with the creation of a file containing the last temperature and salinity profiles. Changes to 
the salinity profile are based upon a numerical solution of die following differential equation describing 
die time dependent behavior of the salinity profile in die absence of any fresh water or brine injection into 
the UCZ. _ 

Eq. (19) 

[dz 


ds m i . d f. 

dt ~ R(z) * dz[ 


•!’•¥]} 


where: 


R(z) = 


the ratio of surface area at depth (z) to the surface area at the top 
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of the NCZ. This term only has significance in small ponds with 
sloping side walls. 

K - the temperature dependent solutal diffusivity of NaCl which 

’ " t*? 0 ' mo,ecular dM *“ k>n of N,a “ 

presence of a salinity gradient. 

s = is the Soret coefficient reflecting the secondary salt transport 

T process that occurs in the presence of a temperature gradient. 

The above etpiatkm has been stapled by considering that th e Sore. ^ “^ NaCl ^io^of 
the order oflxlO» 1* to 3x10- 1* and that for most 

will be relatively small. This was evidenced in analyses performed at Ono State ** « 

o^ed disport differed from the calculated salt transport, using Equation 19 without the Sr term, 

by less than 10% (8]. 

The second stapling adjustment mad. to Equation .9 

diffusivity coefficient. K s , with variations msaluny, “ncennatioos^ toedlfA^ of K ^ 

of ks ^pf^r^r'a r*. z 

po^miri bom d C ma S^om «te Deponent of bnerior's Office of Saiine 

Water 17]. 


K S (I) = (r +95.567 + 2920.9) -9.23 xlO " 


(m 2 thr) 


Eq. (20) 


Finally, the term containing information on the physical aspects of thepondwrils (Le. * W* 

rSUtion isk The geomuric JjVPjr. 

Z^ives, 'from Tmerpolalion of plotted data, a maximum difference insalinny gradient £ 
approximately 10* . The error introduced by excluding this effect ts partially offs ^ y ^ 

Soret coefficient is not included in this analysis. In light of the above considerations 
equation can be reduced to the following, 

Eq. (21) 


i£. V — 

Si ' Si 1 


[ m'-day J 


where: K s is expressed in the solution as a temperature dependent parameter whose 

functional relationship is given by Equation 20. 

This equation was numerically approximated in a manner similar to that used in foVus^in the 
for ^temperature profile. A nodal network, with internodal distances cott^mg to ^ u*d m foe 
determination of the temperature profile, was constructed with nodes placed at the upper and 
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interfaces of the GRZ. The forward difference equations for the salinity profile nodes are given as 
follows: 


Boundary: 



Ea. (22) 


Internal: 



Eq. (23) 


UCZ Node: 




***" t'vcz 

where: S 

as 

the salinity expressed in terms of concentration (i.e., kg/m 3 ) 

K.. 

= 

the solutal diffiisivity expressed in m 2 /Day at the interface 



temperature. 

K., 

= 

the solutal diffiisivity at the upper node’s temperature (m 2 /Day). 

K* 

= 

the solutal diffiisivity at die lower node’s temperature (m 2 /Day). 

Lucz 

= 

the upperzone thickness. 



Discussion of Results 


Two validation runs were performed to assess the accuracy of the computed temperature profile over a 
span of time. A summer month (i.e., August, 1989) and a winter month (i.e., November, 1989) were 
selected and actual measurements, to the extent available, were collected. 

August 1989 Validation Run: Actual measurements taken at the UTEP research pond were used for 

the daily ambient temperature, heat extraction from the storage zone, daily solar insolation, and weekly 
pond profiles for temperature and salinity. The run period extended from August 2, 1989 through August 
31, 1989 for a total of 30 days. A plot of the calculated temperature profile on the last day is plotted 
against actual temperature readings taken on that same day. As can be seen from die comparative plot. 
Figure 4, a close correlation of computed temperature with actual temperature exists throughout the NCZ. 
Only a slight variance between computed and actual storage zone temperatures exists (about .97 e C or 
1.3% from actual). 

November 1989 Validation Run: In this run an initial temperature profile was constructed from pond 

measurements taken on November 7, 1989. During the run period, which spanned 22 days (from 1 1/7/89 
through 1 1/28/89), the salinity gradient was not sufficiently steep at the upper and lower boundaries and 
an adjustment to the thickness of the primary gradient region was necessary. 
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Figure 4. Comparative plot of temperature profiles for 8/3 1 /89. 


comparative plot of the computed temperature profileon the with measurements 

provided in Figure 5. As with the previous Augu* run the , ^ 5 0 . c (U ., 

aken at foe pond. The temperature variance at foe to foe fact that foe LCZ exhibited 

ipproximately 10.0% of acmtal). This higher variance isp during fois period of time. The 

t pronounced temperature stratification (which is D resented in Figure 6. As previously 

esulting salinity profile, derived from taterfaces in order to 

nentioned, a boundary adjustment occurred at both the upper and lower m, 

satisfy foe boundary stability criterion. 



Figure 5. Comparative plot of temperature profiles for 11/28/89. 
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Figure 6. Transition in salinity profile resulting from shifts in NCZ boundaries. 


CONCLUSIONS 

The double diffusive transfer processes typical of a salinity gradient solar pond was represented by a 
system of coupled finite difference equations to simulate tbe interdependency of the ponds temperature 
and salinity gradients. A number of simplifying assumptions can be readily made to allow for an easily 
generated solar pond numerical model that will provide a reasonable first order prediction of changes in 
temperature and salinity profiles, bulk temperature changes in the ponds storage zone, and changes in the 
thickness of the gradient zone. This modeling approach may be used for predicting attainable storage 
temperatures, gradient zone maintenance schedules, and salinity concentration requirements. It is 
applicable to a variety of conditions for purposes of determining die economic viability of a new solar 
pond and could be employed in the operational planning of currently operating solar ponds. The results 
also provide valuable insight into the salinity profiles which need to be installed in order to maintain a 
stable thermal stratification in an established solar pond (for a given storage zone temperature associated 
with a given application). Further, the time variance calculation of the profiles provides useful 
predictions of the salinity gradient modifications which need to be accomplished during die time evolution 
of the operating pond in order to sustain operation. 
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